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LONG-TERM  GOALS 

The  long  term  goals  of  this  project  are  to  investigate  statistical  models  for  signals  propagating  in 
long-range  underwater  channels  and  to  design  signal  processing  techniques  to  mitigate  signal 
fluctuations  due  to  random  disturbances  such  as  internal  waves. 

OBJECTIVES 

At  long  ranges,  broadband  receptions  consist  of  early  ray-like  arrivals  and  a  finale  that  is  best  described 
in  terms  of  the  low  order  modes.  The  energetic  low  mode  signals  are  more  strongly  affected  by  internal 
wave  scattering  than  the  ray  arrivals.  By  focusing  on  the  low  modes,  this  project  seeks  to  develop  a 
better  understanding  of  internal  wave  effects.  The  first  objective  of  this  project  is  to  derive 
range-dependent  mode  statistics  from  experimental  data  obtained  during  the  SPICE04  and  LOAPEX 
experiments.  The  second  objective  of  this  project  is  to  develop  new  robust  signal  processing  techniques 
to  estimate  the  travel  times  of  the  modes  based  on  the  derived  statistical  channel  model. 

APPROACH 

To  characterize  internal  wave  effects  on  the  modes,  this  project  used  the  extensive  data  sets  of 
low-frequency  receptions  recorded  as  a  part  of  the  North  Pacific  Acoustic  Laboratory  (NPAL)  project. 
Two  specific  experiments  are  particularly  relevant  for  the  current  work.  First,  the  Long  Range  Ocean 
Acoustic  Propagation  Experiment  (LOAPEX)  conducted  in  2004  provided  a  unique  opportunity  to 
measure  low  mode  receptions  at  a  series  of  ranges  from  50  km  to  3200  km.  In  addition  to  LOAPEX,  the 
SPICE04  experiment  included  transmissions  from  a  bottom-mounted  source  at  Kauai  to  a  receiving 
array  at  a  range  of  2400  km.  This  project  analyzed  the  LOAPEX  and  Kauai  receptions  and  has 
compared  the  results  to  Parabolic  Equation  (PE)  simulations.  The  results  of  this  analysis  were  used  to 
develop  random  channel  models  for  the  low  order  modes  and  subsequently  develop  new  signal 
processing  techniques  for  these  modes. 
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WORK  COMPLETED 

Mode  Tomography  and  Signal  Processing :  Perturbations  to  the  mode  excitation  spectrum  due  to 
internal  waves  were  analyzed  using  simulated  data,  and  a  model  based  on  second-order  statistics  was 
constructed.  Matched  subspace  detectors  based  on  this  statistical  model  were  used  to  estimate  the  mode 
travel  times  at  LOAPEX  stations  T50  and  T250.  The  estimated  travel  times  were  used  to  invert  for  the 
range- averaged  sound  speed  profile  (SSP)  along  the  LOAPEX  path.  The  inverted  profiles  were 
consistent  with  point  measurements  made  during  the  experiment.  These  results  were  presented  at  the 
North  Pacific  Acoustic  Laboratory  (NPAL)  workshop  in  August  2009  and  the  ASA  meeting  in  October 
2009  [5].  Two  journal  articles  are  currently  in  preparation. 

Internal  Wave  Simulations :  Internal  wave  effects  across  a  400  km  range  were  modeled  by  using  the 
method  suggested  by  Colosi  and  Brown  [6].  In  order  to  analyze  mode  propagation  through  internal 
waves,  two  types  of  simulations  were  implemented.  First,  coupled  mode  simulations  were  performed  to 
calculate  the  amplitude  and  phase  statistics  of  the  unscattered  mode  signals.  Second,  parabolic  equation 
(PE)  simulations  were  used  to  estimate  range-dependent  statistics  of  the  total  mode  signals  (unscattered 
plus  scattered).  Results  of  this  simulation  study  were  presented  at  the  Acoustical  Society  of  America 
(ASA)  meeting  in  November  2006  [3],  the  NPAL  workshop  in  May  2007,  the  ASA  meeting  in  June 
2007  [4],  the  NPAL  workshop  in  May  2008,  and  the  ASA  meeting  in  June-July  2008  [1]. 

LOAPEX  Analysis'.  The  low  mode  signals  received  during  LOAPEX  were  processed.  Mooring  motion 
data  was  missing  for  some  days  of  the  LOAPEX  experiment.  A  method  was  developed  to  estimate  the 
missing  data,  and  those  corrections  were  applied  to  the  receptions  prior  to  mode  processing.  The 
interpolation  method  developed  for  the  mooring  data  was  presented  at  the  MTS/IEEE  Oceans 
Conference  in  September  2008  [2]. 

Kauai  Analysis :  The  signals  received  during  the  2004  NPAL  experiment  from  the  source  deployed  to 
the  north  of  Kauai  were  processed  and  analyzed.  Statistics  of  the  low  order  modes,  such  as  mean, 
variance,  mean  power,  kurtosis,  and  skewness,  were  estimated.  The  results  of  this  analysis  were 
presented  at  the  NPAL  Workshop  in  April  2006. 


RESULTS 


This  section  discusses  signal  processing  methods  for  estimating  the  mode  travel  times  and  the 
implementation  of  tomographic  inversions  using  those  travel  times.  The  first  part  focuses  on  the  use  of 
matched  subspace  detectors  for  travel  time  estimation.  The  second  part  describes  mode  tomography 
using  the  LOAPEX  data  set. 


Matched  Subspace  Detectors'.  This  section  considers  the  problem  of  mode  travel  time  estimation  for 
acoustic  tomography.  Assuming  adiabatic  propagation,  the  received  time  series  for  mode  m  at  a  range  r 
can  be  written 


ej(cot~kmr) 

am{t )  =  ft  I  S{  ft))  (f)m  ( CO,  Zsrc)  /_  dCO, 

'  if 


where  /3  is  a  constant,  5(a))  is  the  source  spectrum,  </>,„(  ft),  zsrc)  is  the  modeshape  at  the  source  depth, 
and  km  is  the  range-averaged  mode  wavenumber.  A  standard  way  to  estimate  the  travel  time  of  a  known 
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signal  is  to  implement  a  matched  filter  (MF)  and  use  the  peak  of  the  MF  output  to  determine  the  arrival 
time.  There  are  several  ways  to  design  the  MF,  depending  on  what  is  assumed  about  the  waveguide. 

The  simplest  way  to  design  the  MF  for  mode  signals  is  to  assume  that  the  modeshape  at  the  source 
depth  is  constant  across  the  band  of  interest  and  that  dispersion  effects  are  negligible.  In  this  case  the 
received  mode  signal  is  a  scaled  version  of  the  source  signal  s(t)  defined  below: 

s(t)=P  I  S((Q)ejmd(0. 

The  MF  for  a  signal  s(t)  is  the  time-reversed  version  of  the  signal,  i.e.,  /?mf  =  s(—, t).  For  the  remainder 
of  this  discussion,  this  version  of  the  MF  is  called  the  standard  MF. 

Since  the  spectrum  of  the  modeshape  at  the  source  depth  is  not  always  constant,  it  makes  sense  to 
extend  the  MF  design  to  incorporate  the  effects  of  </>,„(  CO,  zsre).  Assuming  that  the  modeshapes  are  only 
a  function  of  the  known  background  (quiescent)  sound  speed  profile  and  that  dispersion  effects  are 
negligible,  the  received  mode  signal  is 

^quiescent (0  =  ft  j  >S'( CO)0m( CO,  Zsrc)^'  d(0. 

The  matched  filter  that  incorporates  the  scaling  due  to  the  quiescent  spectrum  of  the  modeshape  at  the 
source  depth  is  flquiescent  =  ^quiescent ( — t).  The  quiescent  MF  is  most  useful  when  the  mode  spectrum 
varies  substantially  across  the  band  of  the  source,  which  typically  happens  when  the  source  is  located 
near  a  null  in  the  mode. 

The  Matched  Subspace  Detector  (MSD),  proposed  by  Scharf  and  Friedlander  [11],  is  an  extension  to 
the  matched  filtering  concept  that  allows  for  received  signals  to  be  characterized  by  a  known  subspace, 
rather  than  requiring  them  to  be  described  by  a  single  vector.  This  project  applied  MSDs  to  the  problem 
of  detecting  the  arrival  of  mode  signals  affected  by  internal  waves.  Internal  waves  perturb  the  arrival 
time  of  the  modes  (by  perturbing  the  mode  wavenumbers)  and  perturb  how  the  mode  is  excited  by 
changing  the  modeshape  </>m(co,zsrc).  Assuming  that  the  modeshape  perturbations  are  random,  the 
second  order  statistics  of  the  received  mode  signal  can  be  used  to  define  the  subspace  used  in  the  MSD. 
For  this  project  the  eigenvectors  of  the  correlation  matrix  for  the  internal-wave-perturbed  mode  signals 
were  used  to  construct  the  subspace.  Numerical  simulations  for  the  LOAPEX  environment  at  50  km 
range  showed  that  the  two  eigenvectors  associated  with  the  two  highest  eigenvalues  of  the  correlation 
matrix  are  sufficient  to  describe  most  of  the  mode  energy.  ID  and  2D  MSDs  were  constructed  for  all  the 
modes  and  used  to  estimate  the  mode  travel  times.  Figures  1  and  2  show  the  output  of  the  two  matched 
filters  and  the  ID  and  2D  MSDs  for  modes  1  and  10  for  two  different  simulated  realizations  of  the 
internal  wave  field.  For  mode  1  the  outputs  of  the  four  different  types  of  filters  are  basically  identical. 
The  mode  1  spectrum  is  essentially  flat  across  the  source  band,  thus  including  the  0,„(co,zsrc)  term  does 
not  have  much  effect.  Any  of  the  four  filters  can  be  used  to  determine  the  arrival  time  of  mode  1 ;  the 
arrival  time  is  estimated  as  the  time  associated  with  the  peak  of  the  filter  outputs.  For  mode  10,  the 
results  are  significantly  different.  In  this  example,  the  source  depth  is  near  a  null  of  mode  10,  thus  the 
mode  excitation  spectrum  has  a  greater  effect  on  what  the  received  times  series  looks  like.  The  location 
of  the  null  changes  with  different  internal  wave  realizations,  and  the  shape  of  the  received  time  series 
varies  with  each  realization.  The  2D  MSD  filter  captures  this  variability,  and  the  output  of  this  filter 
contains  a  single  peak  located  at  the  correct  arrival  time  for  both  internal  wave  realizations.  The 
standard  and  quiescent  matched  filters  and  the  ID  MSD  filter  work  well  for  one  of  the  two  realizations, 
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Mode  1  at  T50 


Mode  10  at  T50 


Seconds  Seconds 

Figure  1:  Outputs  of  the  standard  MF,  quiescent  MF,  and  the  ID  and  2D  MSD  processors  for  modes 
1  and  10  of  a  simulation  at  50  km.  The  plots  show  the  results  for  one  realization  of  the  internal  wave 
field.  For  mode  1  the  four  filters  produce  similar  results  and  any  of  the  filters  could  be  used  to 
estimate  the  mode  1  arrival  time.  For  mode  10,  the  2D  MSD  and  quiescent  MF  could  be  used  to 
estimate  arrival  time  since  they  each  produce  an  output  with  a  single  peak  at  the  correct  time.  The 
standard  MF  and  ID  MSD  processors  have  outputs  with  multiple  peaks,  making  time  of  arrival 

estimation  substantially  more  difficult. 


Mode  1  at  T50  Mode  1 0  at  T50 


Seconds  Seconds 

Figure  2:  Outputs  of  the  standard  MF,  quiescent  MF,  and  the  ID  and  2D  MSD  processors  for  modes 
1  and  10  of  a  simulation  at  50  km.  The  plots  show  the  results  for  a  different  realization  of  the 
internal  wave  field  than  Figure  1.  For  mode  1  the  four  filters  produce  similar  results  and  any  of  the 
filters  could  be  used  to  estimate  the  mode  1  arrival  time.  For  mode  10,  the  2D  MSD  and  the  standard 
MF  could  be  used  to  estimate  arrival  time  since  they  each  produce  an  output  with  a  single  peak  at  the 
correct  time.  The  quiescent  MF  and  ID  MSD  processors  have  outputs  with  multiple  peaks,  making 

time  of  arrival  estimation  substantially  more  difficult. 
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but  not  for  the  other.  These  plots  demonstrate  the  effectiveness  of  the  2D  MSD  filter  for  estimating  the 
arrival  times  of  modes  in  this  environment. 

Consider  the  results  of  travel  time  estimation  using  the  standard  MF  and  the  ID  and  2D  MSD 
processors  for  the  LOAPEX  data  set.  Figure  3  shows  the  travel  time  estimates  for  station  T50  produced 
using  the  four  filters.  The  plot  compares  the  travel  time  estimates  with  the  travel  times  predicted  by  the 
sound  speed  profiles  measured  at  the  receiving  array  and  the  source  location  (approximately  50  km 
from  the  receiver).  For  modes  6,  7,  9  and  10,  which  have  nulls  in  their  excitation  spectra  at  the  source 
depth,  the  ID  MSD  and  standard  MF  have  significantly  larger  travel  time  variance  than  the  2D  MSD,  as 
expected.  The  travel  time  estimates  of  all  three  methods  are  closer  to  the  predictions  for  the  Seabird 
profile  (measured  at  the  array)  than  to  the  CTD  profile  measured  at  the  source.  Figure  4  shows  the  travel 
time  estimates  for  station  T250.  At  this  range,  there  is  not  much  difference  between  the  estimates 
provided  by  the  three  different  methods. 

At  500  km,  the  number  of  eigenvectors  required  to  represent  most  of  the  energy  in  modes  1  to  10  is 
significantly  greater  than  the  number  required  at  T50.  Simulations  at  500  km  showed  that  MSDs  of 
dimensions  ranging  from  4  to  9  were  required  to  describe  most  of  the  energy  in  modes  1  to  10.  Figure  5 
shows  the  travel  time  estimates  obtained  using  the  MSDs  for  the  first  10  modes  at  500  km.  The  MSD 
estimates  are  compared  with  those  obtained  by  picking  the  peaks  of  the  output  of  a  standard  matched 
filter.  The  left  plot  shows  the  mean  of  the  travel  time  estimates,  and  the  right  plot  shows  the  standard 
deviation  of  the  estimates.  The  MSD  travel  time  estimates  are  similar  to  the  travel  time  predictions  of 
the  background  profile.  The  peak-picking  approach  has  a  high  standard  deviation  on  the  order  of 
40-50  ms.  Most  of  the  travel  time  variance  of  peak  picking  is  due  to  the  significant  amount  of 
internal-wave-induced  multipath.  The  MSD  on  the  other  hand  has  a  much  smaller  travel  time  standard 
deviation  on  the  order  of  10-15  ms.  The  standard  deviation  of  the  MSD  is  primarily  due  to  travel  time 
wander  due  to  internal  waves  variations  and  not  due  to  errors  associated  with  scattering.  Based  on  the 
travel  time  simulation  examples  and  travel  time  estimation  statistics  such  as  those  presented  in  this 
section,  MSDs  are  flexible  enough  to  accommodate  the  different  types  of  internal  wave  scattering  for 
each  mode  and  the  different  ranges. 

Mode  Tomography :  This  section  describes  the  implementation  of  mode  tomography  for  the  LOAPEX 
environment  and  then  discusses  the  results  of  the  inversion.  The  inverse  implemented  for  LOAPEX 
used  the  perturbation  theory  approach  of  Munk  and  Wunsch  [9],  where  perturbations  in  mode  travel 
time  are  inverted  for  perturbations  in  sound  speed.  Based  on  previous  work  [7,  8,  12,  10],  it  is  assumed 
that  the  sound  speed  perturbations  can  be  characterized  by  a  small  set  of  depth-dependent  basis 
functions.  For  LOAPEX  the  basis  consisted  of  a  set  of  empirical  orthogonal  functions  (EOF’s)  derived 
from  environmental  data  recorded  on  the  vertical  line  array  used  in  the  experiment.  The  time-averaged 
SSP  at  the  array  was  used  as  the  reference  profile,  and  the  EOF’s  were  computed  from  an 
eigendecomposition  of  the  correlation  matrix  of  the  recorded  perturbations  about  the  reference  profile. 
The  input  to  the  inversion  was  the  mode  travel  time  perturbations,  obtained  by  subtracting  the  mode 
travel  times  for  the  reference  profile  from  the  mode  travel  times  measured  using  the  methods  described 
in  the  previous  section. 

Figure  6  shows  the  results  of  the  tomographic  inversion  at  50  km  range  (LOAPEX  station  T50), 
obtained  using  the  three  types  of  travel  time  estimates  shown  in  Figure  3.  Each  plot  in  the  figure 
compares  the  estimated  sound  speed  perturbation  with  the  perturbation  associated  with  the  CTD 
measurement  at  the  T50  source  location  and  the  Seabird  measurement  at  the  receiver  location.  Since  the 
2D  MSD  method  produces  travel  time  estimates  with  lower  variance,  it  is  not  surprising  that  the 
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Figure  3:  Travel  time  estimates  at  LOAPEX  station  T50  using  the  standard  MF  and  the  ID  and  2D 
MSD  processors.  The  errorbars  indicate  the  standard  deviation  of  the  estimates.  Predicted  travel 
times  are  shown  for  two  measured  sound  speed  profiles  and  a  background  reference  profile.  The 
prediction  denoted  Seabird  is  the  for  the  profile  measured  at  the  receiving  array  and  the  prediction 
denoted  CTD  T50  is  for  the  profile  measured  at  the  source.  The  travel  time  estimates  are  closer  to  the 
Seabird  profile  than  the  source  profile.  The  2D  MSD  has  a  lower  travel  time  variance  than  the  other 

two  methods. 
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Figure  4:  Travel  time  estimates  at  LOAPEX  station  T250.  Errorbars  indicate  the  standard  deviation 
of  the  estimates.  Predictions  are  shown  for  three  measured  profiles:  the  Seabird  profile  measured  at 
the  receiving  array,  and  the  CTD  profiles  measured  at  stations  T50  and  T250.  The  estimates  for  the 
three  methods  all  lie  close  to  the  travel  times  predicted  by  the  measured  sound  speed  profiles  at  T50 
and  T250.  Similar  to  the  results  for  T50,  the  MSD  travel  times  have  a  lower  variance  than  the 

standard  MF  estimates. 
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Figure  5:  Mean  ( left  plot)  and  standard  deviation  ( right  plot)  of  travel  time  estimates  at  500  km  from 
simulations.  The  MSD  travel  time  estimates  have  a  mean  that  is  close  to  the  travel  time  estimates  of 
the  background  SSP.  The  standard  deviations  of  the  MSD  travel  time  estimates  are  much  smaller 
than  the  standard  deviations  of  the  peak-picking  results. 


inverted  profile  obtained  from  the  2D  MSD  processor  has  lower  variance  than  the  ID  MSD  and 
standard  MF  processors.  Note  that  there  is  a  significant  difference  in  the  sound  speed  at  the  source 
(CTD  profile)  and  the  sound  speed  at  the  receiver  (Seabird  profile).  The  inverted  profile  is  closer  to  the 
receiver  measurement  than  to  the  source  measurement,  suggesting  that  the  recorded  modes  propagated 
primarily  in  an  environment  resembling  that  at  the  receiver.  Figure  7  shows  the  results  of  the  inversion 
at  250  km  range  (station  T250)  obtained  using  the  travel  time  measurements  shown  in  Figure  4.  Similar 
to  T50,  the  2D-MSD  result  at  T250  shows  significantly  lower  variance  than  the  other  two  methods.  At 
250  km  range,  the  inverted  profile  lies  between  the  CTD  measurements  obtained  at  stations  T50  and 
T250,  though  it  is  closer  to  the  T250  measurement.  Figure  8  shows  the  inversion  results  for  a  range  of 
500  km  (station  T500).  Figure  8  compares  the  inverted  result  with  the  Seabird  measurement  at  the 
receiver,  CTD  measurements  at  three  LOAPEX  stations,  and  the  average  of  the  Seabird  profile  and  the 
three  CTD  profiles.  The  inverted  result  at  T500  is  closest  to  the  average  profile,  which  makes  sense 
since  the  modes  effectively  propagate  through  the  range- averaged  environment. 

In  summary,  this  project  has  successfully  implemented  mode  tomography  for  ranges  up  to  500  km.  The 
inverted  sound  speed  profiles  show  good  agreement  with  environmental  measurements.  At  ranges  up  to 
500  km,  the  MSD  approach  for  estimating  travel  times  provides  a  significant  reduction  in  variance  of 
the  estimates. 

IMPACT/APPLICATIONS 

This  research  has  signal  processing  and  scientific  applications.  First,  this  work  has  demonstrated  that 
subspace  processing  methods  can  substantially  improve  detection  and  estimation  in  an  environment 
affected  by  internal  waves.  Second,  this  project  has  demonstrated  the  feasibility  of  mode  tomography  at 
ranges  up  to  500  km.  Since  the  low  modes  are  concentrated  around  the  sound  channel  axis,  mode-based 
tomography  should  enhance  the  sound  speed  resolution  at  axial  depths. 


8 


2D  MSD 


ID  MSD 


Standard  MF 


1600  L 


400  - 


600 


1000 


0 

m/s 


1400 


1600 


400 

/* 

\ 

600 

1  ^ 
\ 

/  I 
i  E 

N  H 

800 

\fj 

cf~r 

1000 

r 

1200 

- 

1400 

----- 

1600 

i 

0 

m/s 


Figure  6:  Inverted  sound  speed  perturbations  at  LOAPEX  station  T50.  The  inverted  sound  speeds  are 
close  to  those  associated  with  the  Seabird  measurement  made  at  the  receiving  array.  The  errorbars 
show  that  the  2D  MSD  processor  has  lower  variance  than  the  other  two  processors. 
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Figure  7:  Inverted  sound  speed  perturbations  at  LOAPEX  station  T250.  The  inversion  results  are 
close  to  the  CTD  measurement  at  T250.  The  errorbars  show  that  the  2D  MSD  processor  leads  to  a 

lower  variance  in  the  inversion  results. 
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Figure  8:  Inverted  sound  speed  perturbations  at  LOAPEX  station  T500  obtained  using  the  standard 
MF  travel  time  estimates  for  modes  1  to  10.  The  mean  of  the  sound  speed  inverse  is  close  to  the 
average  SSP  obtained  by  averaging  the  measured  profiles  across  the  T500  path. 
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RELATED  PROJECTS 


This  work  is  closely  related  to  ONR  Award  N00014-09-1-0675,  which  is  the  grant  funding 
Dr.  Chandrayadula’s  advisor  Kathleen  Wage.  It  is  also  related  to  the  work  being  done  by  the  North 
Pacific  Acoustic  Laboratory  group  coordinated  by  principal  investigators  Peter  Worcester  (Scripps)  and 
James  Mercer  (APL  -  UW).  Many  other  ONR-sponsored  researchers  work  on  projects  related  to  NPAL 
and  participate  in  the  NPAL  workshops. 
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